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Abstract 

Modern numerical techniques employing properties of flux Jacobian matrices are ex- 
tended to general, nonequilibrium flows. Generalizations of the Beam- Warming scheme, 
Steger- Warming and van Leer flux-vector splittings, and Roe’s approximate Riemann 
solver are presented for three-dimensional, time- varying grids. The analysis is based on 
a thermodynamic model that includes the most general thermal and chemical nonequi- 
librium flow of an arbitrary gas. Various special cases are also discussed. 

Introduction 

The design of the next generation of advanced space transportation systems such as 
the National Aerospace Plane (NASP) and the Aeroassisted Orbital Transfer Vehicle 
(AOTV) requires detailed flow computations. The combination of high speed and low 
density result in the departure of air from a perfect gas due to the excitation of internal 
modes, dissociation, and ionization. At sufficiently low altitudes (i.e. sufficiently high 
density), the rate of collisions between particles is high enough so that all these processes 
are at equilibrium with their respective reverse processes. The only modification in the 
formulation of the equations is the replacement of the perfect gas law by a general, 
equilibrium gas law [1]. The corresponding extensions of the numerical algorithms are 
found in Ref. 2. There is a significant altitude (i.e. density) range in which the mean 
free path between collisions is sufficiently small for the continuum approximation to 
be valid, but the collision rate is not large enough to maintain thermal or chemical 
equilibrium. This series of papers is devoted to the numerical computation of general, 
nonequilibrium flows. 

Most modern techniques used in CFD utilize the properties of flux Jacobian matrices 
for the treatment of the inviscid terms in the numerical solution of conservation laws. 
For central difference methods, the Beam- Warming scheme [3] requires the true flux 
Jacobian matrices, and their eigenvalues and eigenvectors are needed for the diagonal 
algorithms [4]. Most upwind methods, such as the Steger- Warming flux-vector splitting 
[5], the van Leer flux- vector splitting [6], and the Roe approximate Riemann solver [7], 
all utilize the properties of the flux Jacobian matrix. Their original derivations relied on 
the algebraic simplicity of the perfect gas law. The purpose of this paper is to extend 
these methods to nonequilibrium flows. Only those aspects of the schemes affected 


by the nonequilibrium assumption are covered. Other investigations [8-10] have been 
limited to the treatment of flux Jacobian matrices for chemical nonequilibrium. 

We first discuss the thermodynamic model as it relates to the formulation of the 
inviscid flux. It is general enough to include any type of nonequilibrium flow in an 
arbitrary gas. The exact flux Jacobian matrices, their eigenvalues, and eigenvectors 
are then presented. This is followed by the generalizations of Steger- Warming and van 
Leer flux-vector splittings, the generalization of the Roe average used in Roe’s approx- 
imate Riemann solver. Finally, we discuss the formulations for several special cases 
including different treatments of thermal nonequilibrium and chemical nonequilibrium. 
The analyses are presented for three-dimensional flow with time-varying grids. 

Thermodynamic Model 

We consider a gas of mixture of chemical species which may include neutral or ionized 
atoms or molecules, or free electrons. For each species, statistical mechanics provides 
the molecular basis for deriving the macroscopic equations of state, which relate the 
internal energy and pressure to the density and possibly several temperatures. To 
form the fundamental variable, the partition function, one first requires the quantized 
energy levels and degeneracies derived from quantum mechanics. These are given in 
principle by the eigenvalues of the Schrodinger equation. In practice, they can only be 
determined approximately under some assumptions, and with the aid of experimental 
data. 

A given energy level can be characterized by different types of degrees of freedom. 
Free electrons only have translational degrees of freedom. Atoms also possess electronic 
and nuclear degrees of freedom. Molecules have additionally vibrational and rotational 
degrees of freedom . A thermodynamic model is based on the splitting of the energy 
level into several microscopic modes, each characterized by one or more types of degrees 
of freedom. The physical conditions determine the extent to which such a splitting is 
valid. The macroscopic properties also require the distribution of the species over the 
energy levels. This is given by the Maxwell- Boltzmann distribution characterized by a 
single temperature if the system is in thermodynamic equilibrium. (Note that we do 
not attempt to deal with quantum statistics here.) The phenomena we are studying in- 
volve thermal and chemical nonequilibrium processes which are not strictly describable 
by equilibrium statistical mechanics. Nevertheless, these phenomena ordinarily involve 
systems in what is called local equilibrium. This means we assume that at any instant, 
in the neighborhood of any point in space, the energy level can be expressed as a sum of 
independent energy modes, where the distribution of species over each mode can be ap- 
proximated by a Maxwell-Boltzmann distribution corresponding to some temperature 
defined for that mode. 
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For most gasdynamic applications, the nuclear state is rarely altered. Thus, we 
shall not consider any contribution to the internal energy from the nuclear degrees 
of freedom in our model. Under most conditions, we can also assume that the gas 
mixture consists of weakly interacting particles. The splitting of each energy level into 
a translational mode and an internal structure mode is then completely appropriate. 
From a macroscopic point of view, this is equivalent to stating that the species internal 
energy per unit mass e s is the sum of the average translational energy of random thermal 
motion, £* ran , and the average energy of the internal structure, e\ nt . Here the subscript 
s denotes any particular species in the mixture, whether it is a neutral or ionized atom 
or molecule, or a free electron. A superscript denotes a macroscopic mode of energy. 
This splitting leads to the thermally perfect gas law for each species. Specifically, the 
translational energy is proportional to the translational temperature T s , and is given 


by 


tran 

* M 



( 1 ) 


where R a = R/M , , R is the universal gas constant, and M„ is the molar mass of species 
s. Also, the species pressure p e is given by 


P> — P,R,T , , (2) 

where the species density p t is the the mass of species s per unit volume. 

At this stage, c* nt consists of the internal energy from electronic excitation for atoms, 
and the combined energy of vibration, rotation and electronic excitation for molecules. 
To a somewhat poorer approximation, e\ nt for molecules can be split into the electronic 
excitation energy and the combined energy of vibration and rotation. In a further 
approximation, the latter can be split into separate vibrational and rotational ener- 
gies. This leads to the rigid-rotator, harmonic-oscillator model and other models with 
anharmonic corrections. Each macroscopic energy mode may be characterized by a sep- 
arate temperature or several modes that are approximately in equilibrium with each 
other could be characterized by a common temperature. One may consult the book 
by Herzberg [11] or other textbooks on statistical mechanics for the calculation of the 
internal energy for each mode. It is important to point out that the splitting of e' a nt 
for molecules into separate macroscopic modes can not be justified at higher tempera- 
tures. Rotation produces a centrifugal force that affects the vibration, while vibration 
changes the moment of inertia that affects the rotation. The centrifugal force and the 
vibrational frequency also vary with the electronic state because of different electronic 
configurations. Thus, the splitting may cause very large errors in calculating the in- 
ternal energy at very high temperatures [12]. Recently, Jaffe [13] has proposed a more 
rigorous equilibrium thermodynamic model for diatomic molecules. In his model, for 
each electronic state, and a given rotational quantum number, one determines an effec- 
tive intermolecular potential from which one calculates the allowed energy levels below 
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the potential barrier. The maximum rotational quantum number for each electronic 
state is that for which the effective intermolecular potential first becomes completely 
repulsive. e‘ n< is obtained by summing over these electronic-rotational-vibrational en- 
ergy levels for all allowed vibrational and rotational quantum numbers, and over all 
known electronic states. We have modified his model for obtaining the energy levels 
and have written a general computer program for computing the internal energy [1,12]. 
The results show excellent agreement with those in the JANAF tables [14], which have 
only been tabulated below 6000 K. Although the code operates at a speed of 170 
Mflops on a CRAY 2 computer, the computation is rather tedious and cumbersome, 
and it may involve summing over 20,000 to more than 50,000 energy levels for each 
diatomic species. Thus, for practical purposes, we also provide vectorizable, linear 
search, cubic spline interpolations. Typically, for interpolations from a data base of 
10 OK intervals, the maximum error is less than 0.001% for all the air species except 
O and 1V + , for which the maximum error is less than 0.01% [12]. Jaffe [13] has also 
proposed a nonequilibrium thermodynamic model in which the electronic-rotational- 
vibrational energy level is split into an electronic mode, a rotational mode, a vibrational 
mode, and an interaction mode. The electronic and vibrational modes are character- 
ized by a vibrational temperature, the rotational mode is characterized by a rotational 
temperature, and the interaction mode can be characterized by either the vibrational 
or rotational temperature. Because of the strong interaction between modes, we be- 
lieve that this energy level should not be split, but should be characterized by a single 
internal mode temperature for each species. 

In addition to enabling us to determine accurate species internal energies, the ther- 
modynamic model is also needed to determine the rates of energy transfer between 
different modes, as well chemical rate constants and transport coefficients. The choice 
of an appropriate model is therefore a difficult problem. In this paper, we do not at- 
tempt to judge which model should be used. Instead we formulate our algorithms as 
generally as possible so that they could include all possible models, including those 
that we question. 

We restrict our analysis to mixtures of chemical species in which all the heavy parti- 
cles (i.e. atoms and molecules) are close enough in mass that their translational energies 
can always be characterized by a single translational temperature T. Since the mass 
of the electron is so much smaller than those of the heavy particles, the translational 
energy of free electrons could be characterized by a separate electron temperature T e . 
When this occurs, we will use the subscript s' to denote any heavy particle, and the 
subscript e to denote the free electron. With this notation we can write 

TV = T (3) 

for all s'. Under our assumption of a mixture of thermally perfect gases, the pressure 
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p is given as 


( 4 ) 


P — 5>*..r + p e R e T e . 

a' 

We note that the thermally perfect gas assumption may break down at sufficiently 
high temperatures due to the large orbits associated with high electronic states, or 
at sufficiently high densities. The assumption is also not strictly true for the charged 
species, due to the long range nature of the electric field. 

Under the assumptions of our model, e , can in general be expressed as 

ta = t t t (T) + e'.(T') + e (5) 

Here is defined to be the sum of e‘ ran and the part of c* nt that can be characterized 
by the translational temperature T. t\ is the part of c* n< that can be characterized by 
the electron temperature T e . e‘ is the remaining part of e* nt not in equilibrium with 
either T or T e . In practice, t\ is characterized by a temperature T l t , or possible several 
such temperatures. Note the distinction between cj and e* ran , and also e* and e* nt . In 
order to obtain simple additive relations for the internal energies of the mixture (see 
Eqs. (8) and (9)), we have included the energy of formation from a set of elemental 
species in the definition of e*. Also, the arbitrary constant in the definition of e, must 
be determined in terms of the constants chosen for the elemental species. For the free 
electron, we obtain simply 


4 = 4 = 0> and 4 = \ R ' T '- (6) 

For the heavy particles, the exact division of the energies in Eq. (5) depends on the 
choice of models. For example, e *, could be the energy of the bound electrons, and e’, 
the energy in the vibrational modes, while the energy in the rotational mode would be 
included in ([, . The three temperature model used by Park [15] and Lee [16] is one such 
example, where they considered the special case in which the vibrational energy of all 
the diatomic species are characterized by the same vibrational temperature. The model 
used by Candler and MacCormack [17] is another special case in which e*. = 0 for all 
species, t\, = 0 for the monatomic species, and c‘, for the diatomic species contains 
the vibrational energy of the ground state only. At higher temperatures, where the 
splitting of e* nt into separate modes is not justified, we believe the appropriate model 
is €*, = 0 and e*, = c*?‘. Another possible model would have all the internal energy 
characterized by T or T e , so that e* = 0. 
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Flux Jacobian Matrices for a Nonequilibrium Gas 


To formulate a complete set of conservation equations for a multi-temperature, multi- 
component gas mixture in a thermo-chemical nonequilibrium environment is a nontriv- 
ial task. It requires knowledge in advanced kinetic theory and quantum mechanics. 
The main difficulty conies from the fact that the gas does not consist of structureless 
particles. Also, as the temperature gets higher, the upper excited electronic states 
become populated and some atoms and molecules become ionized. Additional com- 
plexity arises since the particles now interact with each other due to long range forces. 
This can have an important influence on transport, relaxation, radiation and chemical- 
reaction processes. Many existing theoretical approaches have neglected this effect or 
have limited it to a certain degree of complexity. Further research and development 
is required for more rigorous formulations. Under some condition the ionized gas can 
not. be treated as electrically neutral. This gives rise to induced electric fields and ad- 
ditional source terms in the momentum and energy equations. A rigorous formulation 
would require the addition of Maxwell’s equations to our system. Instead of dealing 
with all these complications, in this paper we just focus on some numerical aspects 
of the convective terms. None of the transport, energy transfer, chemical source or 
radiation terms will be discussed here. Details will be reported in the future as part of 
this series on nonequilibrium flow computations. 

The variables defining the state of the gas motion are the mass averaged mixture 
velocity u, the temperatures T and T e , and for each chemical species the density p, 
and e*, (or T‘< ). In order to formulate the conservation equations, it is convenient to 
introduce the internal energy of species s per unit volume 

( 1 ) 

with analogous definitions for e ^, , e *, and l \, . Summing over all species, one can express 
the internal energy of the mixture per unit volume as 

i = e' + ^ £*1', (8) 

where 

S’ 

and = (M) 

S 8 

If all the p $ are known, T and T e can be obtained from e f and c e , although this will 
involve an iterative process if all the e*, or e * are not linear functions of T or T e , 
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respectively. It is also convenient to introduce the overall density of the mixture p , 
given by 

p='5Lp» (10) 

a 

and the species mass fraction a,, defined as 

_ P» 

<*, = — • 

P 

Using Eqs. (10) and (11), one can also define the part of the internal energy in equi- 
librium with T e per unit mass of the mixture as 


(ii) 




(12a) 


and for each species the nonequilibrium part of the 
mass of the mixture as 


e 


i 

9 * 



species internal energy per unit 


(12b) 


Note the distinction between e\, and t*. . The enthalpy per unit mass of the mixture is 
then given by 


, e + p 

h = (13c) 

P 

— [ f »'(-^) + Rs'T] + a e R e T e + £ e + £\,. (136) 

a 1 a ' 


In order to obtain maximum generality and greater compactness, we employ the 
vector approach of Ref. 18. The set of conservative variables per unit volume U can be 
represented compactly by the algebraic column vector 


(14) 


where m = pu is the momentum of the mixture per unit volume, and e = e + |pu • u is 
the total energy of the mixture per unit volume. Here p s and e' a , contain NS and NI 
elements respectively, where NS is total number of chemical species, and N I is total 
number of energy modes which are not characterized by T and T e in the mixture. Note 
that m is a physical vector, while both e and l e are scalars. 
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The calculation of the flux of U across a surface element plays a central role in the 
numerical solution of conservation laws. Let n be the unit normal vector in a positive 
direction to a cell surface in a finite-volume grid, or a coordinate surface in a finite- 
difference grid. If v n is the normal component of the velocity of a time- varying surface, 
and u n = n • u, we can define the normal relative velocity component u' = u n — v n . 
The set of inviscid normal flux components per unit area F n is given by the algebraic 
column vector 




' PsU' ’ 

P 


m v! + pn 

E 

= 

eu' + pu n 

K 


i\,u' 

E e . 


Vn' 


(15) 


where M, P and E denote the normal flux of mass, momentum and energy. Note that 
the above equations do not include the kinetic energy of the free electrons in the last 
row of U , or the electron pressure p t in the corresponding row of F n . If these terms were 
included, the equation set would have to be augmented by the momentum equation for 
the free electrons. In addition to increasing the number of equations, it would require 
modeling the momentum transfer between the free electrons and the heavy particles. 
But to avoid the additional complications arising from the introduction of Maxwell’s 
equations, an approximate form of the electron momentum equation is normally used to 
obtain the induced electric field [19, 15-17]. This results in the presence of the gradient 
of p e in the source term for the electron energy, as well as the overall momentum and 
total energy equations. By excluding the electron kinetic energy and electron pressure 
from our electron energy equation, we must add an additional term involving p e to its 
source term. The effect of the presence of p e in the source terms on shock-capturing 
capabilities remains to be investigated. 

We can define a flux Jacobian matrix operator A satisfying dF n = AdU, using the 
convention that in forming the product of a matrix element with a column vector 
element, a dot product is implied if each element is either a physical vector (e.g. u or 
n) or a tensor (e.g. un or nu). The differential of Eq. (4) takes the form 

dp = <Tdp s > -I- R e T e dp e + J ^2 / p l >R i <dT + p e R e dT e . (16) 

3 1 3 ' 


The species specific heats are defined by 


and 


= 


K 

dT 

(17a) 

df e 

3 

dT e 

(17b) 
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Using Eq. (17), one can rewrite Eq. (16) in terms of the differentials of internal energy 
per unit volume as 

dp = X,dp $ + itde* + k e dV, 


where the pressure derivatives are defined as 


and 


Pt 1 ■R* 1 Q a » .Rj> 

_ T>.'p» >ct v,, ~ E,- a * ,c v’ 

e _ Pt OLfRg 

~ Y,,p'C e v. - ’ 

X,* = — Ke\, — K e e*,, 

X e = ReT t - K e € e e = R e T e (l - |* e ). 


Eqs. (4) and (20) can be combined to derive the identity 

P = P»X» + ««* + «*c e . 


(18) 


(19a) 

(196) 

(20a) 

(206) 


( 21 ) 


Introducing the set of differentials of [/, we obtain the final form 

dp = ^P • u -|- dp$ — ku • dm + «de — k^P de*i + («* - «)de e . 


(22) 


Using Eq. (22), the matrix .4 can then be written as 


A = 


S gj*u ocgitfi 0 

( *u • u + Xr) n — w„u un-fcnu + u'I «n 

( I U • U + Xr - H)Un fln - ««nU + U* -KU n ( K e - k)u u 


0 0 
— Kn (« e — «)n 


-£ e u n 


£ ,' n 

e*n 


0 

0 


6 t i r iu' 

0 


0 


u 


(23) 

where H = h + | u • u is the total enthalpy per unit mass of the mixture, I is the identity 
tensor, and 6 sr , 6,< r < are Kronecker deltas. One can easily verify the homogeneity 
property 

AU = F n , (24) 

which is the direct consequence of the assumption of a mixture of thermally perfect 
gases leading to the form of Eq. (4). 

A has the three distinct eigenvalues 


Aj = u', A 2 = v! + c, and A 3 = u — c, 


(25) 
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where c is the frozen speed of sound given by 


— a aX* + + ( K * — K y! e \' 

8 a* 

(26a) 

— ( K + 1 ) — . 

(266) 


P 


The second form is obtained using Eq. (21). Note that the number of distinct eigenval- 
ues is unaffected by the number of species and the number of nonequilibrium energy 
modes in the system. A] is associated with those conservative variables whose flux is 
purely convective. These are the three types of nonequilibrium variables p 4 , e*, and c e , 
and the tangential component of m. Since Aj is a repeated eigenvalue, if a diagonal 
algorithm is applied, the inversion associated with this eigenvalue can be performed 
simultaneously. This demonstrates one of the advantages of the diagonal algorithms, 
since one only needs to solve three different scalar equations in each direction whether 
the gas is perfect, or in thermodynamic equilibrium, or in thermo-chemical nonequilib- 
rium. 

In order to construct the four types of linearly independent eigenvectors associated 
with Aj, we span the plane normal to n by an arbitrary set of two basis vectors b,, 
and the set of reciprocal basis vector b J , satisfying b, • b- 7 = 6j, where 6\ is the 
Kronecker delta. It follows that b; • n = b J • n = 0. One can then obtain the similarity 
transformation A = RAR ~ i , with the diagonal eigenvalue matrix A defined as 




A 


A = 






(27) 


The right eigenvector matrix R can then be written as 


R = 


S,r 

0 

0 

0 

a. 

0‘s 

u 

cbj 

0 

0 

u + cn 

u — cn 

U - ^ 

K 

c(ubi) 

1 

1-*1 

K 

H -j- cu n 

H — cu n 

0 

0 

6»'r' 

0 



0 

0 

0 

1 

£ e 

£ e J 
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(28) 


while the left eigenvector matrix R 1 can be written as 


R - 1 = 



’C 2 6 ar - 0,(fu • U+ Xr) 

a«KU 

— a,/c 

Q t K 

— a,(/c e - /c) 


— c(b jf • u) 

cb ; 

0 

0 

0 

1 

_£ v(f U ‘ U + Xr) 

e\'K u 

-e\,K 

C 2 S a > r' + e‘,/c 

~C\,{K e - K) 

c 2 

-£ e (fu.U + Xr) 

£*ku 

—e e K 

e*K 

c 2 — £ e (« e — k) 


|(fu-U + Xr-CU n ) 

i(cn-/cu) 

2* 

-\k 

IK"*) 


. |(fU-U + Xr + CU n ) 

-|(cn + «u) 

I* 

-1« 

!(«*-«) - 
(29) 


In actual computations, since the eigenvector matrices are used to operate on algebraic 
vectors or matrices, one probably never needs to form them. By inspecting the structure 
of J? -1 , it is more efficient to express it in a different manner. Eq. (22) can be expressed 
compactly as dp = P T dU in terms of the row vector 


= [ *U • U + \, —KU 


K 


— K 


K e — K 


Then R 1 can be partitioned as 




c6 sr 

0 

0 

0 

0 - 

0 


— • u 


0 

0 

0 

-£ e 

1 

P T + ~ 

C 

0 

0 

0 

0 

0 

0 

cS a f 
0 

0 

c 

] 

2 


2 

n 

2 

0 

0 

0 

1 

2 J 


it,, 

L 2 

n 

2 

0 

0 

0 . 


(30) 


(31) 


When carrying out the implicit inversions in a diagonal algorithm, or forming the 
limiters in some TVD schemes, one needs to perform the operations R~*V and RW, 
where V and W are some algebraic column vectors of appropriate structure. Let 


V = 


v 2 

n 

Va,‘ 


L v 5 J 


(32) 


be a five-component algebraic column vector with the structure of U or F n . In terms 
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of the scalar quantities 


and 


Vl 

=E 

g 

Vl., 



(33a) 

V4 

= £ 

Va,', 



(335) 


a 

P T V 




V 6 

C 2 






_ 1 
c 2 

_ a 

Vl. + K ^ 

~(u • u)vj - u • v 2 + v 3 - v^j + (« e - k)v s 

,(33c) 

v 7 

1, 

= -(n • v 2 - 

«nVi), 


(33d) 


one can calculate R 1 V efficiently as 


R~'V = 


- a,V 6 

1 [b J ' • v 2 - (b> *u)t>i] 
V 4 ,' - £ 1 ,>v 6 
vs - e*v 6 
|(r 6 + t’7 ) 

|(t>6 - v^) 


(34) 


Note that each additional species s involves only an additional add in forming rj, 
an additional add and multiply in forming v$, and an additional add and multiply in 
forming the first component of R~ 1 V. Similarly, each additional species s' that contains 
i involves only an additional add in form Va and an additional add and multiply in 
forming the third component of R~* V. Thus the addition of new species to the system 
involves very few additional operations. 

Let 


W = 


Wl» 

U>2 

V>3»' 

W 4 
W 5 

L Wfi J 


(35) 


be any six-component algebraic column vector with structure of R 1 V . In terms of the 
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scalar quantities 


and 


one can calculate RW efficiently as 


RW = 


w i = 

3 

(36a) 

W 3 = Y^ W 3s', 

g 1 

(366) 

U>7 = U>5 + w 6, 

(36c) 

W 6 = W S — W6, 

(36d) 

+ a s w 7 

1 

1 


(ti>i + )u + c(w } 2 bj + u> 8 n) 

X»U>U + c{w } 2 U • b j + U> 8 U n ) + Hw 7 + U >3 + (1 - ^)w 4 

U’ 3j > + £*,U)7 
W 4 -|- t e W-i 


37) 

Again, each additional species s involves an additional add in forming , an additional 
add and multiply in forming the first component of RW, and an additional add and 
multiply in forming the third component of RW. Similarly, each additional species s' 
that contains e’, involves an additional add in forming 103 , and an additional add and 
multiply in forming the fourth component of RW . Thus, also in this case the addi- 
tional of new species to the system involves very few additional operations. On some 
occasions, a matrix of the form RD (/(A/)) J ? -1 may be needed, where I? is a diagonal 
matrix whose elements are function of the eigenvalues. One should never form R and 
i ? _1 individually and then perform the matrix multiplications. The aforementioned 
techniques should also be applied in this case. 

In diagonal algorithms with operator splitting, one needs to form the operator R^Ri, 
where the subscript l refers to a surface with unit normal n j and a set of basis vectors 
bjj, while the subscript k refers to a surface with unit normal n* and a set of reciprocal 
basis vectors b*. From Eqs. (28) and (31) one obtains the matrix 


Rl'Ri = 


[6.r 

0 

0 

0 

0 

0 1 

0 

K • bii 

0 

0 

M 

~K • 

0 

0 

6,'r' 

0 

0 

0 

0 

0 

0 

1 

0 

0 

0 

fn* b u 

0 

0 

|(1 + n* • nj) 

|(1 - n* -nj) 

. 0 

-§n* bu 

0 

0 

|(1 - n k -m) 

§(1 + n * • n*). 


(38) 


Note that the entries involve only the grid geometry, and are independent of the physical 
flow variables. Therefore the presence of chemical or thermal nonequilibrium does not 
require any additional work in operating with the matrix R^ 1 Ri on an algebraic vector. 
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Generalized Steger- Warming Flux- Vector Splitting 


If |«'| > c, it follows from Eq. (25) that the eigenvalues of A are all of one sign, thus 
allowing upwind differencing to be simply implemented. If |tt' | < c, the eigenvalues of 
A are of mixed sign. In flux-vector splitting methods, the flux F n is written as 


F„ = F++F-, (39) 

so that the split-flux Jacobian operators A* satisfying dF± = A ± dU have the property 
that Aj( J 4 ± )^0 for all i. (In practice, one can permit some A t ( J 4 ± )^0 if they are 
sufficiently small in magnitude.) 

Since the homogeneity property Eq. (24) is valid, the two ways of obtaining Steger- 
Warming flux-vector splitting for a perfect gas can be easily extended here. One way, 
the so called “eigenvalue splitting” [2, 20, 21], is to expand F n as 


3 

F n = y, F -> 
1=1 


(40) 


where each F n , is associated with the distinct eigenvalues A*. One can readily obtain 
the expressions 


where 


F ni = 


*1(7-1) 

7 


Ps 


Ps 

P u 

A 2 

p ( u ± cn) 

e " 

?s' 

, F 2 = — 

’ n 3 27 

p{H ± cu n ) 

i e 


l e J 



= K + 1 . 


(41) 


(42a) 

(42b) 


The generalization of Steger- Warming flux-vector splitting is obtained by letting F± 
be the sum of those F n i associated with A>0. For — c < u' < 0 we therefore have 

F+ = F n 2 and F~ = F nl 4 - F n3 , (43) 

while for 0 < u' < c 

F+ = F ni + F n2 and F~ = F n3 . (44) 
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The other way, the so called “plus-minus” splitting, follows the original derivation of 
Steger and Warming [5] for the perfect gas. One can express the eigenvalue A ; as 

+ A i , (45a) 


where 


_ K + \^i\ _ K ~ |^t| 

’ ~ 2 ’ A * “ 2 


(45b) 


Using the homogeneity property and applying the similarity transformation, one ob- 
tains 


FZ = R\+R-'V and F~ = RA~R~ 1 U, (46) 


where A + and A~ have the diagonal elements A + and A“ respectively. Both approaches 
give identical results. While one can not prove rigorously for the nonequilibrium case 
that all the eigenvalues of A* have the appropriate sign, this has been shown [20] to be 
true for equilibrium flow of a diatomic gas with vibrational excitation. It is reasonable to 
assume that this condition will be approximately satisfied for practical nonequilibrium 
flows. Furthermore, the numerical stability depends also on other conditions. 

Due to the fact that A ± are discontinuous when A* changes sign, “glitches” or oscilla- 
tions have been observed at sonic points, stagnation points, and shocks in many perfect 
gas applications. There are many techniques which have been reported in the literature 
to overcome such phenomena. For example, Steger and Warming [5] suggested adding 
small blending terms to A f. This yields a smooth sonic transition. Thomas et. al. [22] 
indicated that by using the MUSCL type differencing [23] one can obtain a smooth 
sonic transition and stable shock structure. In addition, Buning [21] and Ying [24] 
utilized a proper “transition operator” so that the sonic glitches and shock oscillations 
can be avoided. All these techniques can be easily applied to the nonequilibrium flow 
cases. 


Generalized van Leer Flux-Vector Splitting 


van Leer [6] constructed a different type of flux-vector splitting for a perfect gas in 
terms of polynomials of u' whose A* are continuous and have one zero eigenvalue. This 
results in a sharper capt ure of transonic shocks. The generalization for a nonequilibrium 
flow is similar to that derived previously for an equilibrium flow [2,20]. For |u'| < c, 
the continuity conditions necessitate a factor (u' ± c) 2 in the formulas for F*. Using 
the fact that the critical values of u' are ±c, one can easily demonstrate that for any 
f(u') that is either an even or odd function of u', 

/+(«') = ±/-(-«') « f(u')=±f(-u'). (47) 
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As a function of u' , all components of F(U) can be represented by at most the cubic 
polynomials 

f(u') = ao + a\u' + a. 2 U 12 + a 3 u' 3 . (48) 

The lowest order polynomials for f±(u') satisfying Eqs. (39) and (47), as well as the 
continuity conditions at |u'| = c can then be expressed as 

/*(«') = (a* + afu'), (49) 

where 

Oj = ± (aj — c 2 a 3 ) (50a) 

c 

and af — 2 ca 3 ± (a 2 — -r). (506) 

c l 


Using Eqs. (48)-(50), one readily obtains the relations 


M, ± = ±^(u’±c) 2 , 

Ef = ±^u'±c)\ 

E'* = ±?-(u' ± cf , 
Ac 


and 


t ± , pW ± c? 


P 1 = i- 


4c 


1 


u (u' ^ 2c)n 

7 


(51) 

(52) 

(53) 

(54) 


One can similarly obtain a two term expansion for E ± . But this could not reduce 
to van Leer's solution for a perfect gas where A * has one zero eigenvalue. However, 
the function ±{v' 4- c) 2 (u' — c) 2 satisfies the continuity conditions at u' = ±c, and 
such a term can be added to E ± without affecting E(U). For a perfect gas, there is 
only one species, and 7 is a constant. The zero eigenvalue condition results if E ± = 
/(JV/i.P* , n,r„). Guided by the form of Eq. (54), we express the additive function in 
such a way that E ± can be written as 


E ± = ± 


K±c) 2 r {(7 — 1)1/ ± 2c> 2 
4c P 2(7 2 — 1) 


+ e 


7-1 


pu 


1 2 


u 1 T 2c) + rapin' =F c) 2 


(55) 


where V’ is an arbitrary parameter that is independent of the arbitrary constant in the 
definition of e. For a perfect gas, van Leer’s condition results if -0 = 0. It is shown in 
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Ref. 20 that V’ = 0 is also a reasonable assumption for an equilibrium real gas. For the 
general nonequilibrium case, it is also simplest to take = 0. Again, one can not prove 
rigorously that all the eigenvalues of A* have the appropriate sign, but it is reasonable 
that this is approximately true, and one of the eigenvalues will be close to zero. 

Generalized Roe's Approximate Riemann Solver 

Among the various approximate Riemann solvers, the most common one uses the Roe 
average because of its simplicity and its ability to satisfy the jump conditions across 
discontinuities exactly. The derivation in Ref. 8 for a perfect gas employed parameter 
vectors. To obtain a generalization for a nonequilibrium flow, we use a different, more 
direct approach, which is similar to that derived priviously for equilibrium flow [25]. 

In approximate Riemann solvers based on local linearization, the flux at a surface 
separating two states t/j,and Ur is based on the eigenvalues and eigenvectors of some 
average A. The optimum choice for A is one satisfying 

A F n = AAU, (56) 

where A(-) = (• )r — (*)j> This choice of A captures discontinuities exactly. One way 
of obtaining A is to seek an average state U , such that 

A = A(U). (57) 

The notation U implies only those variables that appear explicitly in Eq. (57). Such 
a state is known as a Roe-averaged state, and was derived by Roe for a perfect gas. 
The generalization to a nonequilibrium flow is obtained by substituting Eqs. (14), (15), 
(23), and (57) into Eq. (56), with n and r n fixed at the surface. In general, ul, u r, 
and n are arbitrary, independent vectors. We can therefore expand u as 

u = a L u L + crUr + a„n. (58) 

The independence of \ii and uj? and n also implies the independence of the dot products 
ui -ul, Ul-ur, ur-ur, \iL-n and uj?-n. After substituting Eq. (58) into the momentum 
component of Eq. (56), and equating coefficients of independent quantities, one readily 
establishes that 


y/pL 

(59a) 

aL ~ Ve l + y/ei' 

_ , _ _ y/pR 

&R — 1 &L j ? 

yfPL + y/PR 

(595) 

and o n = 0. 

(59c) 
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Therefore Eq. (58) becomes 

u = o £ ui, -1-aflUfl. (60) 

This is the identical relation derived by Roe for a perfect gas. It satisfies all the terms 
involving the velocity. The remaining terms in the equation result in the new condition 

+ k'Ac* = Ap, (61) 

3 

where Ac* is obtained from the jump of the conservative variables using Eq. (8). This 
is just the discrete form of Eq. (18), averaged between the two states. This condition 
is automatically satisfied for a perfect gas. Substituting Eq. (60) into the first, fourth 
and fifth components of Eq. (56), one obtains the relations 

a, = aia,i + aRa,R, 

e\, = a L e\. L + a R e \, R , 
and £ e = &i,£ e l + a R£ e R- 

The third component of Eq. (56) results in the additional relation 

H = a L HL + clrHr, (65) 

which is also true for a perfect gas. From the definition of H , Eqs. (60) and (65) can 
be combined to define a Roe-averaged specific enthalpy as 

h = aihi + ORfiR + -atOflAu • Au. (66) 

Note that h could lie outside the range given by hi and Hr if the magnitude of Au is 
sufficiently large. The Roe-averaged sound speed is given by Eq. (26a) as 

c 2 = ^ X. + Hh + (#c* - k)!* - k £*,. (67) 

S 3* 


(62) 

(63) 

(64) 


For an arbitrary nonequilibrium flow, Eq. (61) provide only one relation for the 
variables \ s , k, and n e . We thus have the paradoxical situation that not only does a 
Roe-averaged exist for a nonequilibrium flow, its precise value is not uniquely defined. 
For the special case in which L and R are precisely those that satisfy the jump conditions 
across a discontinuity, Eqs. (59) through (67) are consistent with the exact Riemann 
solver, even though x 4 , k, and K e are not uniquely defined. For a gas dynamic shock 
wave, which necessitates the conditions Aq 4 = Ae*, = At e = 0, one obtains [25] 


h = 


a L + a R 


( 68 ) 
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and 


c 2 


A p 
A ? 


(69) 


The values of h and c 2 as given by Eqs. (68) and (69) will in general not be consistent 
with Eqs. (13b) and (26a). 


There are two approaches to obtain unique values of \ a , 7c, and K e . In one approach, 
average temperatures are defined which are then substituted into Eqs. (19) and (20). 
The other approach is based on thermodynamic states L and R, and does not attempt 
to define an average thermodynamic state. In the first approach one can also distinguish 
between those that employ h giving by Eq. (66), and those that do not. Since h depends 
on the velocity u l and ujj, the resulting values of , 7c, and k e will also depend on 
these velocities. 


In the first approach, one defines e* in terms of a, and T e using Eq. (12a) and solves 
for T e using the Roe conditions Eqs. (62) and (64). tt e and \ e are then obtained from 
Eqs. (19b) and (20b). If we do not use h , Eqs. (19a) and (20a) can be considered as 
parametric definitions of a curve in the \»'~ K space. A unique solution for \ t > and k 
is obtained if the curve has one intersection with the hyperplane defined by Eq. (61). 
Unfortunately, there may be cases where there is no intersection, or possibly multiple 
intersections. On the other hand, one can also obtain T from h using Eq. (13b). The 
resulting unique values of and k derived from Eqs. (19a) and (20a) will not, in 
general, satisfy Eq. (61). 

The second approach is an extension of the method described in Ref. 25. Eqs. (9a) 
and (19a) implicitly define the function /c(p 4 <, e t ). One similarly defines K e (p s ,i e ) from 
Eqs. (9b) and (19b), and the functions x*(p 4 , e*,( e ) from Eqs. (9) and (20). The straight 
line path between states L and R in the space are defined parametrically by 


p»{r) = p t L + tA/j 4 , (70a) 

( t (T) = € i L + rAe t , (70 b) 

V(t)=V l + tL?, (70 c ) 

where the parameter r is chosen so that tj, = 0 and tr = 1. The integrated averages 

X»= [ X* [/M T M‘( T M e ( T )] dr, (71a) 

t/0 

7c = / K [Mt),«*(t)] dT, (716) 

Jo 

and K e — [ K e [/>«(t),c'(t)] dr (71c) 

Jo 


satisfy Eq. (61) exactly. 
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The use of Eqs. (71) is not practical, so that some approximate quadratures are 
required. The resulting values of K i and k' will also not satisfy Eq. (61) exactly. 
Let x»i « and « e define either approximations to Eq. (71) or the values derived from 
h in the first approach. One then requires values of x,i and k' satisfying Eq. (61) 
which are closest to x,, k and k'. This can be formulated geometrically as projecting 
the point x«? * and k' onto the hyperplane defined by Eq. (61). But in order for the 
Roe- averaged quantities to be independent of the arbitrary constant in the definition of 
c a , one must first recast the problem so that geometric relationships will not be affected 
by the choice of this constant. This can be accomplished if one first divides Eq. (61) 
by k. The orientation of the hyperplane in the space defined by the variables x 4 /Jc, 
1 /«, and K e /k is now uniquely defined by states L and R. A further scale factor d with 
the dimension of x , must be introduced, since the x* are not dimensionless. The final 
relations can be written as 


_ _ Dx» + g 2 A p, 6p 
Xt ~ D - Ap6p ’ 

__ Dk 

D — A p Sp ’ 

_ Dk' + A l'8p 
K ' ~ D-ApSp ’ 

where 

bp = A p — ^ Xs Ap, — k Ac* — k' l' 
8 

and D = ^(g Ap,) 2 -I- (Ap) 2 + (Ac*) 2 . 

a 

A natural choice for the scale factor d is 


(72a) 

(726) 

(72c) 

(73) 

(74) 


S = ?. (75) 

The simplest approximation to Eqs. (71) is the trapezoidal rule 

p: _ (-)l + (•)« 

U 2 

although other approximations, such as middle point rule or Simpson’s rule, can be 
used. Note that the magnitude of the error 6p using Eq. (76) is a function of the 
nonconvexity of the variables x»i K^and k' between the two states. On the other hand, 
the magnitude of bp derived from h is also a function of Au. 

An important quantity in the approximate Riemann solver is the column vector 
R~ 1 AU. Its components are the jumps in the characteristic variables. While the 


(76) 
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above relations are all that are required to construct a Riemann solver using Roe’s lin- 
earization, an additional algebraic simplicity can be achieved by expressing differences 
in conservative variables in terms of differences in primitive variables. If one formally 
defines 

P ~ y/ P lPRi (77) 

the expression can be written more simply as 


R-'AU 


Ap, - a t Ap/c 2 
~phj • Au/c 
Ac’, - e\,Ap/c 2 
Ac* — e'Ap/c 2 
\ (A p/c 2 + pn • Au/c) 

\ (A p/c 2 - pn • Au/c) _ 


(78) 


Special Cases 

The results presented so far are for the most general case of nonequilibrium flow. 
There are a number of special cases, which can be derived from the general case by 
deleting one or several equations, or modifying some of the terms in the equations. 
These will be detailed below. 

We first consider the cases relating to the treatment of thermal nonequilibrium. One 
such case occurs when c e = 0, but l\, ^ 0. In practice, this would only occur if there 
were no ionization. Then s and r would equal s' and r' in all the equations, and p e , 
( e s , , e e and their differentials would be set equal to zero in all the equations in which 
they appear. In the algebraic vectors and matrices, the fourth row would be absent 
in Eqs. (27), (29), (31), (34), (35), (38), and (78), while the fifth row would be absent 
in Eqs. (14), (15), (23), (28), (32), (37), and (41). Similarly, one would eliminate the 
fourth column in Eqs. (27), (28), and (38), and the fifth column in Eqs. (23), (29), 
(30), and (31). Note that it is possible to have nonequilibrium resulting from internal 
structure alone even if there is only one chemical species present . 

The other special case of thermal nonequilibrium occurs when c*, = 0 and ( e ^ 0. In 
practice, this would only occur with non- negligible ionization, so that p e ^ 0. Under 
these assumpt ions, the nonequilibrium part of c*? 1 would all be characterized by T e . For 
this case, all e*, and their differentials would be set equal to zero in all the equations in 
which they occur. In the algebraic vectors and matrices, the third row would be absent 
in Eqs. (27), (29), (31), (34), (35), (38), and (78), while the fourth row would be absent 
in Eqs. (14), (15), (23), (28), (32), (37) and (41). Similarly, one would eliminate the 
third column in Eqs. (27), (28), and (38), and the fourth column in Eqs. (23), (29), 
(30), and (31). 
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Both of the above cases can be combined if all the internal energy l s is at equilibrium 
with the translational temperature T, and we only have chemical nonequilibrium. Then, 
e' j( , c e , and their differentials would be absent in all the equations. Since any free 
electrons would have the temperature T, the distinction between the heavy particles 
and the free electrons is no longer valid. Thus s' would be replaced by s in all the 
relations, and terms involving T e would be absent in Eqs. (4), (13b), and (16). Both 
sets of rows and columns listed in the previous two paragraphs would be deleted from 
their respective algebraic vectors and matrices. 

We next consider the cases relating to the treatment of chemical nonequilibrium. 
There are two cases where p e is not present as a conservative variable in the definition 
of U. One occurs when we assume that charge neutrality exists locally at every point. 
In terms of the species ionic valence Z t >, one can express the free electron density as 


Pe — r ^ ] Ps'R, 1 ^ s’ 


(79) 


Eq. (79) is used in evaluating Eqs. (4), (9b), (10), (12a), (13b) and (19b). Eq. (10) can 
be written as 

P ~ fit' Pa' j (80) 

s' 


where 


13s' 


= 1 + 


R 3 ' Z s > 
Re ' 


(81) 


From Eq. (11) it follows that 

= 1 - (82) 


The major change in the equations is to replace s and r by s' and r' in Eqs. (14), (15), 
(18), and every equation starting with Eq. (21). It is also necessary to replace Eq. (20a) 
by 

XV = RAT + Z,r e ) - Kt\, - « e (c*, + ^ Rs'Z.'Te ). (83) 

For most engineering purpose, one can take /? a . ss 1, since R a ’Z s >/R c <1. If this 
small term is not neglected, then the factor /? 4 < must appear in a number of terms. 
This results from expressing dp and du 1 in terms of the differentials of the conservative 
variables. The term u-u in Eqs. (22) and (30) must be multiplied by /?,». The following 
terms must be also multiplied by (3 r i : all the terms in the first column of A in Eq. (23) 
except the u' in the first row and the Xr' in the second and third row; the terms 
involving u in the first column of R in Eq. (28); all the terms involving u and u n in 
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the first column of R 1 in Eq. (29) and the second matrix of Eq. (31). In addition, 
Eqs. (33a) and (36a) must be replaced by 

V, = ^ 

s' 

and 

s' 

The other case in which p t is not a conservative variable is in the absence of ionization, 
where p e = 0. This can be considered as a specialization of the previous case to Z t > = 0. 
It follows from Eq. (81) that /3 t > = 1 rigorously, so that individual terms need not be 
modified. The main difference is that we must also take i e = 0, with the accompanying 
changes described earlier. 

There are also cases which involve only one conservation equation of mass. The 
subscripts s, r, s', and r' are then absent. One such case, described earlier, is thermal 
nonequilibrium due to the internal structure for a gas consisting of a single species. 
A more important case is that of a mixture in thermodynamic equilibrium, for which 
e 1 = e e = 0. Strictly speaking, the flow can not be in truely thermodynamic equilibrium 
unless the flow is uniform. Nevertheless, under the assumption of local equilibrium, all 
the thermodynamic variables can be determined from the process of maximization of 
the irreversible entropy. Thus one can replace the species density equations by a single 
global density equation. All the results for this case have been presented by the authors 
in an earlier paper [2]. However, one can also deduce the results from the general forms 
presented here if one formally replaces a by 1. All the relations are still valid, except 
that the equation of state is given by 


(84) 

(85) 


p = pOm)> 

and Eqs. (19) and (20) are replaced by 

**(!),■ 

The other change is that Eqs. (26b) and (42b) are no longer valid. The equation of state 
and the pressure derivatives, Eqs.(86)-(88), can be obtained from Ref. 1. Similarly, for 
the perfect gas case, one needs to replace Eqs. (86)-(88) by 


p = KC, 

(89) 

k = const., 

(90) 

X = 0- 

(91) 


( 86 ) 

(87) 

( 88 ) 
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As one can see, the main difference in the equations for nonequilibrium, equilibrium 
and perfect gas flows is in the expressions for the pressure and its derivatives. The 
difference for one-dimensional, two-dimensional, or three-dimensional flows is just the 
number of components in the vectors u, m, n, b,, and b J . Thus one can easily build a 
universal computer code which does nonequilibrium, equilibrium and perfect gas flow 
computations in either one, two, or three dimensions. 


Conclusions 

A thermodynamic model has been established for the most general thermal and chem- 
ical nonequilibrium flow of an arbitrary gas. Extensions of modern CFD techniques to 
general, nonequilibrium flows have been obtained. The results have been presented in 
a form that reduces readily to various special types of nonequilibrium flows, as well as 
equilibrium and perfect gas flows. One can therefore construct a universal code that 
covers all these types of flows in any number of dimensions. 
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